/*********************************************************************
*                                                                    *
* inidat -- provides initial data                                    *
*                                                                    *
* Parameters:                                                        *
* y_a_ptr        -- pointer where to write the initial data          *
* initial_time   -- initial time                                     *
*                                                                    *
* Returns: pointer to field_array where data was writen              *
*                                                                    *
*********************************************************************/

#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
#include "gen_2d.h"
#ifdef SV
#include "java_ser.h"
#include "java_fser.h"
#endif



/* struct field_array *inidat(struct field_array *y_a_ptr) { */

void inidat(struct field_array *y_a_ptr, 
	    struct GRID_PAR *grid_ptr,
	    struct INI_PAR *ini_par_ptr){

/* -------> grid parameters <--------------------------*/


  
    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 

    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;

/*     int grid = (*grid_ptr).grid; */

#ifdef EXCISION
    int exc_n1 = (*grid_ptr).exc_grid_1;
    int exc_n2 = (*grid_ptr).exc_grid_2;
    int exc_n1i = (nf_1-ni_1-1)/2-exc_n1+1; /* take out a -1 for visual compatibility with the rest of code */
    int exc_n2i = (nf_2-ni_2-1)/2-exc_n2+1;
    int exc_n1f = (nf_1-ni_1-1)/2+exc_n1; 
    int exc_n2f = (nf_2-ni_2-1)/2+exc_n2;
#endif

    FLOAT x_i = (*grid_ptr).initial_x;
    FLOAT x_f = (*grid_ptr).final_x;
    FLOAT y_i = (*grid_ptr).initial_y;
    FLOAT y_f = (*grid_ptr).final_y;
    FLOAT swap, sigma_b02, r2;

    int coorx = 1;
    int coory = 2;

/* -------> initial data parameters <------------------*/

 FLOAT twoPIdN1 = 2.*PI*(x_f-x_i)/(FLOAT)(nf_1-ni_1-1);
 FLOAT one_dN1 = (x_f-x_i)/(FLOAT)(nf_1-ni_1-1);


 FLOAT twoPIdN2 = 2.*PI*(y_f-y_i)/(FLOAT)(nf_2-ni_2-1);
 FLOAT one_dN2 = (y_f-y_i)/(FLOAT)(nf_2-ni_2-1);

 
  /* Parameters from main */


  FLOAT a0 = (*ini_par_ptr).a0;                
  FLOAT k_a_10 = (*ini_par_ptr).k_a_10; 
  FLOAT k_a_20 = (*ini_par_ptr).k_a_20;   
  FLOAT shift_a0 = (*ini_par_ptr).shift_a0;          

/* Amplitude sin(k_a1* x + shift_a1) in U1 */

  FLOAT a1 = (*ini_par_ptr).a1;                
  FLOAT k_a_11 = (*ini_par_ptr).k_a_11;  
  FLOAT k_a_21 = (*ini_par_ptr).k_a_21; 
  FLOAT shift_a1 = (*ini_par_ptr).shift_a1;          

/* Amplitude of cos(k_a0* x + shift_a0) in U0 */

  FLOAT c0 = (*ini_par_ptr).c0;                
  FLOAT k_c_10 = (*ini_par_ptr).k_c_10;    
  FLOAT k_c_20 = (*ini_par_ptr).k_c_20; 
  FLOAT shift_c0 = (*ini_par_ptr).shift_c0;          

/* Amplitude of cos(k_c1* x + shift_c0) in U1 */
  FLOAT c1 = (*ini_par_ptr).c1;                
  FLOAT k_c_11 = (*ini_par_ptr).k_c_11; 
  FLOAT k_c_21 = (*ini_par_ptr).k_c_21;
  FLOAT shift_c1 = (*ini_par_ptr).shift_c1;          

/* Amplitude of exp(cos(k_b0*x)^2/sigma_b0) in U0 */
  FLOAT b0 = (*ini_par_ptr).b0;                
  FLOAT sigma_b0 = (*ini_par_ptr).sigma_b0;
  FLOAT c0_1 = (*ini_par_ptr).c0_1;
  FLOAT c0_2 = (*ini_par_ptr).c0_2;

/* Amplitude of exp(cos(k_b1*x)^2/sigma_b1) in U1 */
  FLOAT b1 = (*ini_par_ptr).b1;                
  FLOAT sigma_b1 = (*ini_par_ptr).sigma_b1;
  FLOAT c1_1 = (*ini_par_ptr).c1_1;
  FLOAT c1_2 = (*ini_par_ptr).c1_2;

/* Global wave motion */

  FLOAT v1 = (*ini_par_ptr).v1;
  FLOAT v2 = (*ini_par_ptr).v2;

  int initial_data_type = (*ini_par_ptr).initial_data_type;


  /* --------------------------------------------------------------------------------------- */

#ifdef DEBUG_INIDAT
  printf("Begining of function inidat\n");
#endif

#ifdef EXCISION
{
if ((FLOAT)exc_n1*one_dN1 > 0.38*m) 
printf("excision region too large, R_exc_1 = %f \n", (FLOAT)exc_n1*one_dN1);

if ((FLOAT)exc_n2*one_dN2 > 0.38*m)
printf("excision region too large, R_exc_2 = %f \n", (FLOAT)exc_n2*one_dN2);

}
#endif

  /* --------------------------------------------------------------------------------------- */

  /* The cross-like coordinates */

switch ((*grid_ptr).grid)
{
 case 3:
  swap = x_i - (x_f-x_i);
  x_f = x_f - (x_f-x_i);
  x_i = swap;
  break;
 case 4:
  swap = y_i + (y_f-y_i);
  y_f = y_f + (y_f-y_i);
  y_i = swap;
  break;
 case 5:
  swap = y_i - (y_f-y_i);
  y_f = y_f - (y_f-y_i);
  y_i = swap;
  break;
 default:
  swap = x_i + (x_f-x_i)*(*grid_ptr).grid;
  x_f = x_f + (x_f-x_i)*(*grid_ptr).grid;
  x_i = swap;
  break;
}


/*---------> values for different fields <-------------*/


 

/* struct field_array y; */

/* first the time */

/* y.a.time = initial_time; */

(*y_a_ptr).time = (*grid_ptr).initial_time;



  
  /*----------- Initial data for the scalar field ----------- */

  switch (initial_data_type){
  case 0: // give data using cross-like coordinates with fourth order polynomial



      {register int g_ind1, g_ind2;
   FLOAT x=0.;
   FLOAT y=0.;
   FLOAT r=0.;

  for (g_ind1 = ni_1; g_ind1 < nf_1; ++g_ind1) {
    for (g_ind2 = ni_2; g_ind2 < nf_2; ++g_ind2) {
	  x =  x_i+(FLOAT)(g_ind1)*one_dN1;
	  y =  y_i+(FLOAT)(g_ind2)*one_dN2;

/*     (*y_a_ptr).u[PHI][IND_2(g_ind1,g_ind2)] = 	 */
/*        (a0*sin(  (FLOAT)g_ind1*twoPIdN1*k_a_10  */
/* 		 + (FLOAT)g_ind2*twoPIdN2*k_a_20  */
/* 		 + shift_a0) */
/*      + c0*cos(  (FLOAT)g_ind1*twoPIdN1*k_c_10  */
/* 		+ (FLOAT)g_ind2*twoPIdN2*k_c_20 */
/* 		+ shift_c0) */
/*        )*b0*exp(-sigma_b0*sigma_b0*( */
/* 		  (x-c0_1)*(x-c0_1) */
/* 		  + (y-c0_2)*(y-c0_2) */
/* 		  ) */
/* 		); */

	  r2 = (x-c0_1)*(x-c0_1) + (y-c0_2)*(y-c0_2);
	  sigma_b02 = sigma_b0*sigma_b0;
	 if (r2 <= sigma_b02){
	     (*y_a_ptr).u[PHI][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),4)/pow(sigma_b02,4); 
	 }
	 else
	     {(*y_a_ptr).u[PHI][IND_2(g_ind1,g_ind2)] = 0.0;}
    }
  }
      }

    

    break;

  case 1: // l=2,m=2 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0+(FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0+(FLOAT)(grid_ind2)*one_dN2; 
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((1.+x*x)/(1.0+x*x+y*y))*
		(2.0*(1.0/(1.0 + x*x))-1.0);
	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
		
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
			
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((1.+x*x)/(1.0+x*x+y*y))*
		(2.0/(1.0 + x*x)-1.0);
	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((1.+x*x)/(1.0+x*x+y*y))*
		(2.0*(x*x/(1.0 + x*x))-1.0);
	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
			
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((1.+x*x)/(1.0+x*x+y*y))*
		(2.0*(x*x/(1.0 + x*x))-1.0);
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((x*x+y*y)/(1.0+x*x+y*y))*
		(2.0*(y*y/(y*y + x*x))-1.0);
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =   -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =   -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
	
	    if(x*x+y*y>-0.000001){			    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((x*x+y*y)/(1.0+x*x+y*y))*
		(2.0*(y*y/(y*y + x*x))-1.0);
            }
	    else{
              (*y_a_ptr).u[PHI]GRD_IND = c0*(-0.5+1.5/(1.0+x*x+y*y))*
                (2.0*(1.0)-1.0);
	    }     	
	    }}}
	break; 
      }
    
    break;
    
  case 2: // give data using cross-like coordinates with 8th order polynomial 
      {register int g_ind1, g_ind2;
   FLOAT x=0.;
   FLOAT y=0.;
   FLOAT r=0.;

  for (g_ind1 = ni_1; g_ind1 < nf_1; ++g_ind1) {
    for (g_ind2 = ni_2; g_ind2 < nf_2; ++g_ind2) {
	  x =  x_i+(FLOAT)(g_ind1)*one_dN1;
	  y =  y_i+(FLOAT)(g_ind2)*one_dN2;

	  r2 = (x-c0_1)*(x-c0_1) + (y-c0_2)*(y-c0_2);
	  sigma_b02 = sigma_b0*sigma_b0;
	 if (r2 <= sigma_b02){
	     (*y_a_ptr).u[PHI][IND_2(g_ind1,g_ind2)] = 
		 b0*pow((r2-sigma_b02),8)/pow(sigma_b02,8); 
	 }
	 else
	     {(*y_a_ptr).u[PHI][IND_2(g_ind1,g_ind2)] = 0.0;}
    }
  }
      }
    break;

/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
/* CALSE l=2,m-1 multipole */
/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
  case 3: // l=2,m=1 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*(y/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 
			
	      (*y_a_ptr).u[PHI]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2, grid_ind3;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*(-y*x/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
	      (*y_a_ptr).u[PHI]GRD_IND = c0*((y*x)/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	
	      (*y_a_ptr).u[PHI]GRD_IND = c0*(-y/(1.0+x*x+y*y));
	    }}}
	break; 
      }
    
    break;


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/

/* CALSE l=4 multipole */

  case 4: // l=4 multipole
    {register int grid_ind1, grid_ind2;
		    
    for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
      for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
	  (*y_a_ptr).u[PHI]GRD_IND = 0.0;
	}}}
    break;

  default: //unknown initial data type
    printf("scalar field, unknown data type\n %d", initial_data_type);
    exit(1);
    break;
  }


  /*--------- data for the x derivative of u[0] --------------*/


DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[PHI],(struct field *)&(*y_a_ptr).u[PHI_X],&coorx);

  /*---------- data for the y derivative of u[0] ------------- */
DERIV(grid_ptr,(struct field *)&(*y_a_ptr).u[PHI],(struct field *)&(*y_a_ptr).u[PHI_Y],&coory);



  /*----------- data for the time derivative of u[0] ---------------- */
  switch (initial_data_type){/* initial_data_type switch */

    case 0: {// give data using cross-like coordinates
      {register int grid_ind1, grid_ind2;
      FLOAT x=0.;
      FLOAT y=0.;

      for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {

	    x =  x_i+(FLOAT)(grid_ind1)*one_dN1;
	    y =  y_i+(FLOAT)(grid_ind2)*one_dN2; 

	    (*y_a_ptr).u[PHI_T]GRD_IND = 	
	      (a1*sin(  (FLOAT)grid_ind1*twoPIdN1*k_a_11 
			+ (FLOAT)grid_ind2*twoPIdN2*k_a_21
			+ shift_a1)
	       + c1*cos(  (FLOAT)grid_ind1*twoPIdN1*k_c_11 
			  + (FLOAT)grid_ind2*twoPIdN2*k_c_21
			  + shift_c1)
	       )*b1*exp(-sigma_b1*sigma_b1*(
					    (x-c1_1)*(x-c1_1)
					    + (y-c1_2)*(y-c1_2)
					    ))
	      + 
	      v1*(*y_a_ptr).u[PHI_X]GRD_IND
	      + 
	      v2*(*y_a_ptr).u[PHI_Y]GRD_IND;
	  }}
      }}
      break;
  
    case 1: {// l=2 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		(-0.5 + 1.5/(1.0 + x*x + y*y));
	    }}}
	break;
				
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		(-0.5 + 1.5/(1.0 + x*x + y*y));
	    }}}
	break;

      default: // this if for patches 0 to 3
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 

	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		(-0.5 + 1.5*y*y/(1.0 + x*x + y*y));
	    }}}
	break; 
      }
    }

      break;
    

    case 2: {// give data using cross-like coordinates
      {register int grid_ind1, grid_ind2;
      FLOAT x=0.;
      FLOAT y=0.;
      FLOAT z=0.;
      FLOAT r = 0.;

      for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
	    x =  x_i+(FLOAT)(grid_ind1)*one_dN1;
	    y =  y_i+(FLOAT)(grid_ind2)*one_dN2; 
	    r = (x-c1_1)*(x-c1_1) + (y-c1_2)*(y-c1_2);

	    if(r >= sigma_b0){
	      (*y_a_ptr).u[PHI_T]GRD_IND = 0.0;
	    }

	    else{
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		  pow(r-sigma_b1,8)/pow(sigma_b1,8)    
		+ 
		v1*(*y_a_ptr).u[PHI_X]GRD_IND
		+ 
		v2*(*y_a_ptr).u[PHI_Y]GRD_IND;
	    }
	  }}}
    }
      break;
    


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
/* CALSE l=2,m-1 multipole */
/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/
  case 3: // l=2,m=1 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 0:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(y/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 2:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 1:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y*x/(1.0+x*x+y*y));
	    }}}
	break;
				
      case 3:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
			
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*((y*x)/(1.0+x*x+y*y));
	    }}}
	break;
		
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));
	    }}}

	break;
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
	
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*(-y/(1.0+x*x+y*y));
	    }}}
	break; 
      }
    
    break;


/* -------------------------------------------------------------*/
/* -------------------------------------------------------------*/


    case 4: {// l=4 multipole
      switch ((*grid_ptr).grid){ // patches 4 and 5 are the ones that have the poles
      case 4:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		0.125*(3.0 - 30.0/(1.0 + x*x + y*y) + 35.0/(1.0 + x*x + y*y)/(1.0 + x*x + y*y));
	    }}}
	break;
				
      case 5:
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
				    
	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		0.125*(3.0 - 30.0/(1.0 + x*x + y*y) + 35.0/(1.0 + x*x + y*y)/(1.0 + x*x + y*y)); 
	    }}}

	break;

      default: // this if for patches 0 to 3
	{register int grid_ind1, grid_ind2;
	FLOAT x=0.;
	FLOAT y=0.;
			
	for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1) {
	  for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2) {
				    
	      x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2; 

	      (*y_a_ptr).u[PHI_T]GRD_IND = c1*
		0.125*(3.0 - 30.0*y*y/(1.0 + x*x + y*y) + 35.0*y*y*y*y/(1.0 + x*x + y*y)/(1.0 + x*x + y*y));
	    }}}
	break; 
      }
    }

      break;
 
    default: // unknown initial data type
      printf("unknown initial data type = %d \n", initial_data_type);
      exit(1);
      break;

    } 




  printf("<LI>Inidat finished </br></LI>\n"); 
}




